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ABSTRACT 

Ultracompact and hypercompact H II regions appear when a star with a mass 
larger than about 15 solar masses starts to ionize its own environment. Recent obser- 
vations of time variability in these objects are one of the pieces of evidence that suggest 
that at least some of them harbor stars that are still accreting from an infalling neutral 
accretion flow that becomes ionized in its innermost part. We present an analysis of 
the properties of the H II regions formed in the 3D radiation-hydrodynamic simula- 
tions presented by Peters et al. as a function of time. Flickering of the H II regions 
is a natural outcome of this model. The radio-continuum fluxes of the simulated H II 
regions, as well as their flux and size variations are in agreement with the available 
observations. From the simulations, we estimate that a small but non-negligible frac- 
tion 10 %) of observed H II regions should have detectable flux variations (larger 
than 10 %) on timescales of 10 years, with positive variations being more likely to 
happen than negative variations. A novel result of these simulations is that negative 
flux changes do happen, in contrast to the simple expectation of ever growing H II 
regions. We also explore the temporal correlations between properties that are directly 
observed (flux and size) and other quantities like density and ionization rates. 
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1 INTRODUCTION 

The most massive stars in the Galaxy, 0-type stars 
with mas ses Mt, > 20 Mq , emit copious amounts of UV 
photons (|Vacca et al.l 1 19961 ) that ionize part of the dense 
gas from which they form. The resulting H II regions are 
visible via their free-fre e continuum and recombina tion line 
radia tion (|Mezger fc H enderson 1967; Wood & C hurchwelll 
1 19891 ). H II regions span orders of magnitude in size, from gi- 
ant (D ~ 100 pc) bubbles, to "ultracompact" (UC) and "hy- 
percompact" (HC) H II regions, loosely defined as those with 
sizes of ~ 0.1 pc and ~ 0.01 pc (or less), respectively (sec the 
reviews by [Churchwell2002. : ,Kurta.2005. : .Hoare et al.,,2007. ') . 
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UC and HC H II regions are the most deeply embedded, and 
so are best observed at radio wavelengths. 

Large, rarefied H II regions expand without interruption 
within the surrounding medium due to the high pressure 
contr ast between the ionized and neutral phases (|Spitzeil 
|l978j ). This simple model was extrapolated to the ever 
smaller objects recognized later and is widely used to in- 
terpret observations of UC and HC H II regions. Common 
assumptions about these objects are: i) They are steadily 
expanding within their surrounding medium at the sound 
speed of the ionized gas, ~ 10 km s~^. ii) The ioniz- 
ing star(s) is already formed, i.e., accretion to the massive 
star(s) powering the H II region has stopped. 

However, evidence has accumulated that suggests a re- 
vision of these assumptions: 

(i) The "hot molecular cores" embedding UC and 
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HC H II regions are often ro t ating and infalling (e.g. , 
Ketdll990l: ICesaroni et allllQQSl : iBeltran et aLlbOOd I2OI0I : 



Klaassen et al.ll2009l l. sometimes from parsec scales all the 



way to the immediate surroundings of the ionized region 
l|Galvan-Madrid et al.ir2009l : iBaobab Liu et al.ll2010l ). 

(ii) Infall of gas at velocities of a few km s~^ directly 
toward the ionized center has also been observed in UC 
and HC H II region s (Zhang fc H0II1997I : iBehran et al]l2006l : 
iGalvan-Madrid et al. 2009 ). 

(iii) The inner ionized gas has been resolved in a few 
cases, and it also s hows accr etion dynamics (outflow, infall, 
and rotation, 



_e^ 



.Keto fc W ood 2006; Scwilo ct al. 2008; 
"al.l l2009l ). 

(where the flux goes 
H II regions is ~ 



iGalvan-Madrid et 

(iv) The spectral index a (wliere the Hux goes as 
Sv oc v") of some UC and HC H II regions is ~ 1 
from cm to mm wavelengths, indicating density grad ients 
and /or dumpiness inside the ionized gas iFranco et al.ii2000l: 
Ignace fc Churchwelll |2004 iKeto et al. I I2OO8I : lAvalos et al.l 



2009^ ■ 



(v) A few UC and HC H II regions h ave been shown to 
have variatio ns on timescales o f years j Acord et al.l 199^; 
Franco-Hernandez fc Rodrfgued | 2004l : van der Tak et al] 
20051 : ICalvan-Madrid et all |2008| ). These variations in- 
dicate t hat UC an d HC H II regions sometimes 
expand (lAcord et al.| Il998l) . and sometimes contract 
iGalvan-Madrid et al.l I2OO8I ). Some other ionized regions 
around massive protostars h ave been shown to remain ap- 
proximately constant in flux ()Goddi et al.ll201ll '). 

All these observations strongly suggest that UC and 
HC H II regions are not homogeneous spheres of gas freely 
expanding into a quiescent medium, but rather that these 
small H II regions are intimately related to the accretion 
processes forming the massive stars. Simple analytic models 
show that the observed H II region can be either the ionized, 
inner part of the inflowing accretion flow (Kcto 2002, 2003) 
or th e ionized photoev aporative o utflow ()Hollenbach et al.i 
1 1994 ) fed by accretion l|Ketdl200'i1 ). 

Numerical simulations of the formation and expansion 
of H II regions in accretion flows around massive stars 
have only recently become possible with three-dimensional 
radiation-hydrodynamics. Studies that simulate the expan- 
sion of H II regions have focus sed on larg er-scale effects 



on the parental rn o lecular cloud ijPale et a l 2005, ' 2007a 3 
Peters et al.) I2OO8I : ICritschneder et all 12009 : Arthur et all 



201 j ). However, in none of those simulations was the ionizing 



radiation produced by the massive stars dynamically form- 
ing through gravitatio nal collapse in the m olecular cloud. 
Recent simulations by IPeters et al.l ()2010al ) (hereafter Pa- 
per I) include a more realistic treatment of the formation 
of the star cluster. Radio-continuum images generated from 
the output of the simulations of Paper I show time variations 
in the morphology and flux from the H II regions produced 
by the massive stars in formation. These changes are the 
result of the complex interaction of the massive filaments 
of neutral gas infalling to the central stars with the ionized 
regions produced by some of them. A statistical analysis of 
the H II region morphologies (jPeters et al ] |2010bl . Paper II) 
consistently reproduces the relatively high fraction of spher- 
ical and unresolved regions found in observational surveys. 
Thus, the non-monotonic expansion of H II regions, or flick- 
ering, appears able to resolve the excess number of observed 



UC and HC H II regions with respect to the expectation 
if th ey expand uninterrupted (the so-called lifetime prob- 
lem, IWood fc Churchwelll I1989I ). Furthermore, the ionizing 
radiation is unable to stop protostellar growth when accre- 
tion is strong enough. Instead, accretion is stopped by the 
fragmentation of the gravitationally unstable accretion flow 
in a process we call "fragmentation-induced s tarvation", a 
theoretical discussion of which can be found in lPeters et al.l 
(2010c) (Paper HI). 

In this paper we present a more detailed analysis of 
the flux variability in the simulations presented in Paper 
I. In §2 we describe the set of numerical simulations and 
the methods of analysis. In §3 we present our results. §4 
discusses the implications of our findings. In §5 we present 
our conclusions. 



2 METHODS 

2.1 The Numerical Simulations 

Our study uses the highest resolution simulations of 
those presented in Paper I. The simulations use a modi- 
fied v ersion of the adaptive-mesh code FLASH (|Frvxell et al.l 
l2000l ). including self-gravity and radiation feedback. They 
include for the first time a self-consistent treatment of gas 
heating by both ionizing and non-ionizing radiation. We re- 
fer the reader to Paper I for further details of the numerical 
methods. The initial conditions are a cloud mass of 1000 
Mq with an initial temperature of 30 K. The initial density 
distribution is a fiat inner region of 0.5 pc radius surrounded 
by a region with a decreasing density cx r~^'^ . The density 
in the homogeneous volume is 1.27 x 10"^'^ g cm~^. The sim- 
ulation box has a length of 3.89 pc. The size and mass of the 
cloud are in agreement with those of star-formation regions 
that are able to produce at least one star with hh > 20 Mq 
(e.g., Galvan-Madrid ct al. 2009). 

Run A (as labeled in Paper I) has a maximum cell res- 
olution of 98 AU and only t he first collapsed obje ct is fol- 
lowed as a sink particle (see iFederrath et al.ll20ld ). In this 
run, the formation of additional stars (the sink particles) 
is suppressed using a density-dependent temperature fioor 
(see Paper I for details). On the other hand, in run B ad- 
ditional collapse events are permitted and a star cluster is 
formed, each star being represented by a sink particle. The 
maximum resolution in Run B is also 98 AU. 



2.2 Data Sets 

For the entire time span of Run A (single sink) and Run 
B (multiple sink) , radio-continuum maps at a wavelength of 
2 cm were generated from the simulation output every ^ 300 
yr by integrating the radiative transf er equation for free-free 
radia tion while n eglecting scatter i ng ( Gordon fc Sorochenkd 
I2OO2I ) . Following iMac Low et al.l (|l99ll ) , each intensity map 
was then convolved to a circular Gaussian beam with half- 
power beam width HPBW — 0.14" (assuming a source dis- 
tance of 2.65 kpc). A noise level of 10"'^ Jy was added to 
each image. Further details are given in Paper I. These maps 
were used to explore the behavior of the free-free continuum 
from the H II region over the entire time evolution of the 
simulations. For Run B, sometimes the H II regions overlap 
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both physically in space and/or in appearance in the line 
of sight. Therefore, the presented time analysis refers to the 
entire star cluster unless otherwise specified. 

To compare more directly to available observations, 
which span at most a couple of decades in time, each of 
Run A and B were re-run in four time intervals for which 
a flux change was observed in the radio-continuum images 
mentioned above. Data dumps and radio-continuum maps 
were generated at every simulation time step (~ 10 yr). The 
analysis performed in the low time-resolution maps was also 
done in these high time-resolution data. For Run B, the in- 
tervals for the high time-resolution data were also selected 
such that the H II region powered by the most massive star 
is reasonably isolated from fainter H II regions ionized by 
neighboring sink particles, both in real space and in the 
synthetic maps. 



3 RESULTS 

3.1 Variable HII Regions 

Real UC and HC H II regions, as well as those that 
result from the s imulations pr esented here, are far from the 
ideal cases (e.g.. lSpitzeilll978l ). However, it is instructive to 
discuss the limiting ideal cases to show that their variability 
is a natural consequence of their large to moderate optical 
depth. 

The flux density Si, of an ionization-boundecQ H II re- 
gion is 



(1) 



where k is the Boltzmann constant, c is the speed of light, 
V is the frequency, and the brightness temperature Tb is 
integrated over the angular area Q of the H II region. 

In the limit of very low free-free optical depths {rg <C 1), 
Tb along a line of sight / goes as 



Tb ocTr^'^^t/"^ ' / n'^dl, 



(2) 



where Te and n are the electron temperature and density 
respectively. 

Combining equations (1) and (2) we have that for a 
given f and constant T^: 



S{ts < 1) oc n'^R^ oc iV, 



(3) 



where the last proportionality comes from the Stromgren 
relation N (x n {N is the ionizing-photon rate and R is 
the 'radius' of the H II region). Equation (3) shows that in 
the optically-thin limit the H II -region flux only depends on 
A'^. For time intervals of a few x 10 yr the mass and ionizing 
flux of an accreting protostar remain almost constant, and 
so does the flux of an associated optically- thin H II region. 
On the other hand, for very high optical depths {rg ^ 



^ For H II regions that are embedded in their parental cloud, 
the ionization-bounded approximation is better than the density- 
bounded approxi mation. The analysis here presented can b e 
derived from, e.g., iMezger &: H enderson (1967); Spitzcrl l ll97Sll ; 
iRvbicki fc LightmanI l ll979l) : and lKetol ((200^). 
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Figure 1. 2-cm flux (5'2cm, flUed black squares) and stellar mass 
(M*, gray circles) as a function of time (t) for Run A. Although 
the long-term trend of the H II region is to increase in flux, it 
constantly flickers during its evolution. The fluxes at the different 
projections (X-axis top, Y-axis middle, Z-axis bottom), though not 
the same, follow each other closely. 



1) Tb — Te- For a given frequency and constant Te, equation 
(1) becomes: 

S{ts » 1) (X oc A2/^n-^/^ (4) 

therefore, the flux of an optically thick H II region is pro- 
portional to its area, and both flux and area decrease with 
density. 

This analysis is valid for time intervals larger than the 
recombination timcscale (~ 1 month for n ~ 10^ cm~^, see 
e.g., .Osterbroch, 19891 ) and as long as the growth of N is 
negligible. 

However, the H II regions in the simulations are clumpy 
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Figure 2. 2-cm flux (52cmi filled black squares) of the H 11 region 
formed by the most massive star and its stellar mass (M*, gray 
circles) as a function of time (t) for Run B. Although the long- 
term trend of the H II region is to increase in flux, it constantly 
flickers during its evolution. Only the Z-axis projection is used 
because from this viewing angle, perpendicular to the plane of 
the accretion flow, the brightest H II region can be distinguished 
from other H II regions at all times. 



and have subregions of high and low free- free optical depth. 
Their flux during the accretion stage (while they flicker) is 
dominated by the denser, optically thicker {rg > 1) subre- 
gions, so their behaviour is closer to cq. (4) than to eq. (3). 
The on-line version of this paper contains a movie of t// 
for Run B as viewed from the Z-axis (line of sight perpen- 
dicular to the plane of the accretion flow). The dumpiness 
and intermediate-to-large optical depth of these H II regions 
arc also the reasons behind their rising spectral indices up 
to relatively large frequencies {v > 100 GHz) without a sig- 
nificant contribution from dust emission (see the analytical 
discussions of Ignace & Churchwell 2004 and Keto et al. 
2008, for an analysis of these simulations see Paper II). As 
for the variability, the largo optical depths cause the size and 
flux of the simulated H II regions to be well correlated with 
each other, and anticorrelated with the density of the cen- 
tral ionized gas (see Section 3.6). The neutral accretion flow 
in which the ionizing sources are embedded is filamentary 
and prone to gravitational instability (further discussion is 
in Section 3 of Paper III, see also Paper II). The changes in 
the density of the H II regions are a consequence of their pas- 
sage through density enhancements in the quickly evolving 
accretion flow. 



3.2 Global Temporal Evolution 

Figure 1 shows the global temporal evolution of the H 11 
region in Run A (single sink particle). The 2-cm flux (yS'2cm) 
observed from orthogonal directions and the mass of the 
ionizing star (Af*) arc plotted against time. The global tem- 
poral trend of the H II region is to expand and become 
brighter. However, fast temporal variations are seen at all 
the stages of the evolution. The fluxes in the projections 
along the three different cartesian axes follow each other 
closely. For the rest of the analysis, the Z-axis projection, a 
line of sight perpendicular to the plane of the accretion flow, 
is used. 

The H II region is always faint (yS2cm < 1 Jy at the 
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Figure 3. 2-cm flux variations (A52cin) for Run A {top panel) 
and Run B {bottom) panel. 



assumed distance of 2.65 kpc) for M* < 25 A4q- Past this 
point, the H II region is brighter than 1 Jy 86 % of the time 
(Fig. 1). 

A similar analysis of the H 11 region around the most 
massive star in Run B (multiple sinks) is shown in Fig. 2. 
Only the Z-axis projection, i.e., a line of sight perpendicu- 
lar to plane of the accretion flow is used, since only from 
this viewing angle the brightest H II region is well sepa- 
rated from other H 11 regions at all times (the flux movies 
are presented in Paper I). The extra fragmentation in Run B 
translates into a weaker accretion flow and a lower-mass ion- 
izing star as compared to that of Run A, a process referred 
to as 'fragmentation-induced starvation' (a theoretical dis- 
cussion of this process in presented in Paper HI). Therefore, 
the brightest H II region in Run B (Fig. 2) is weaker than 
the H II region in Run A (Fig. 1). Figures 1 and 2 do not 
show times later than t — 100 kyr in order to facilitate their 
comparison. Both runs continue past this time, but in Run 
B accretion onto the most massive star stops &t t — 109 kyr, 
while in Run A the artificial suppression of the fragmenta- 
tion leads to an unrealistically large mass for the ionizing 
star at later times (see Paper 1). 

H 11 regions are highly variable both in Run A and Run 
B. However, since the suppression of fragmentation in Run 
A produces a larger accretion flow and a most massive ioniz- 
ing star, Run A presents larger flux variations than Run B. 
Figure 3 shows the flux changes over the evolution of the H 
II regions in both runs. A comparison of the fractional vari- 
ations of the H II regions shows that, though more similar 
between runs, they are still larger in Run A (Figure 4) . For 
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Figure 4. Fractional flux variations (AS2cm/S'2cm,o) for Run A 
(top panel) and Run B {bottom) panel. 



consecutive data points, in Run A positive flux variations 
are 56 % of the events and the flux increment is +74 % on 
average, while negative changes are 44 % of the events and 
have an average magnitude of —27 %. Similarly, for Run B, 
positive changes (52 % of the events) have an average mag- 
nitude of +42 %, while the average flux decrement (48 % of 
the events) is —18 %. 

3.3 Comparison to Surveys 

We present a comparis on with the ionized-ga s surve ys of 
UC and HC H II re gions bv lWood fc Churchwelll (|l989l ) and 
iKurtz et al. Figure 5 shows normalized histograms 

of 2-cm luminosity (^cmd^) obtained for both runs using 
the time steps previously show n in figures 1 to 4 and the 
CO- a dded obs erved samples from lWood fc Churchwelll (Il989l ) 
and iKurtz et al. ( 1994), taking the 81 sources for which they 
report a 2-cm flux and a distance. Simulation and observa- 
tions roughly agree, but neither Run A nor B can reproduce 
the high-luminosity end of the observed distribution: 20 % 
of the observed sources have S2c-md? > 50 Jy kpc'^, while 
only 1 % of the H II regions in the simulation steps of Run 
A and in no step in Run B have luminosities above this 
threshold. These bright UC H II regions likely correspond 

^ The relation of the surrounding molecular gas to the ionized 
gas has not been explored in detail for many of their sources, 
which makes difficult to assess whether they are candidates to 
harbor accreting protostars. 
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Figure 5. Histograms of 2-cm luminosities {S2cmd?) for the 
global temporal evolution of the simulated H II regions (dashed 
lines) and the co-added samples of Wood & Churchwell (1989) 
and Kurtz et al. (1994) [solid lines). The top and bottom frames 
correspond to Run A and Run B respectively. 



to stages in which accretion to the ionizing protostar(s) is 
completely shut off, therefore they do not correspond to the 
analyzed H II regions from the simulation in which the pro- 
tostars are still accreting. In the range 20 < 5'2cmrf'^ < 50 
Jy kpc^. Run A matches better than Run B the observed 
luminosities. This may be because Run B does not produce 
any star with a mass > 30 AIq . Simulations identical to 
Run B produce higher-ma ss stars in the presence of magnetic 
fields ijPeters et al.ll201lh and may also produce more mas- 
sive stars in the purely radiation-hydrodynamical case with 
higher-mass initial clumps. We plan to perform in the future 
a study on the robustness of our results for different initial 
conditions. For the smallest luminosities. Run B matches 
better than Run A the observed 55 % of H II regions that 
have 5'2cmd^ < 5 Jy kpc^, but over-estimates (35 %) the ob- 
served fracti o n (14 % ) of H II regions with 5 < S2cmd^ < 10. 
iPeters et al] (|2010bh presented statistics of the morpholo- 
gies of the H II regions in Run A and Run B and found that 
Run B agrees better with observed surveys. Although Run 
A can be interpreted as a mode of isolated massive star for- 
mation, the treatment of fragmentation is more realistic in 
Run B (see Paper I). Moreover, mos t high-mass stars form 
in clusters jZinnecker fc Yorkell2007l ). 
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Figure 6. Probabilities for flux increments Icirger than 10 (top 
panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the long-term evolution of the H II region in Run 
A. The error bars indicate the Icr statistical uncertainty from the 
number of counts in each bin 1000-yr wide. 
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Figure 7. Probabilities for flux increments larger than 10 (top 
panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the long-term evolution of the H II regions in Run 
B. The error bars indicate the Icr statistical uncertainty from the 
number of counts in each bin 1000-yr wide. 



3.4 Long-Term Variation Probabilities 

We address the question of the flux variability expected 
from the simulations by calculating the probability of vari- 
ations larger than a given threshold as a function of time 
difference between steps. The low time-resolution data has 
the advantage of spanning the entire runs, but is not useful 
to predict the expected variations on timescales shorter than 
10^ yr. Wc use the high time-resolution data sets to make 
an estimate of the flux variations over shorter timescales, 
but we caution that the analyzed time intervals may not be 
representative of the entire simulation. We use this approach 



because re-running the entire simulations to produce data 
at ^ 10 yr resolution is not feasible. 

Figures 6 and 7 show the probabilities oiflux increments 
larger than a given threshold for time difl^erences between 1 
and 60 kyr for Run A and Run B respectively. The top panels 
correspond to flux increments larger than 10 %, the middle 
panels correspond to 50 %, and the bottom panels to 90 %. 
On average, H II regions tend to expand, making a given 
flux increment to be more likely to happen for larger time 
intervals than for shorter ones. 

Figures 8 and 9 show the probabilities of flux decre- 
ments larger than a given threshold (in modulus) for Run A 
and Run B respectively. At At > 30 kyr the negative-change 
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Figure 8. Probabilities for flux decrements larger than 10 (top 
panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the long-term evolution of the H II region in Run 
A. The error bars indicate the Icr statistical uncertainty from the 
number of counts in each bin 1000-yr wide. 
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Figure 9. Probabilities for flux decrements larger than 10 (top 
panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the long-term evolution of the H II regions in Run 
B. The error bars indicate the Icr statistical uncertainty from the 
number of counts in each bin 1000-yr wide. 



probabilities decrease and reach ~ at about At > 40 kyr 
for any threshold. This is caused by the eventual growth 
of the H II regions in spite of the flickering. There is some 
indication that the probabilities for negative changes also 
decrease at timescales shorter than 1 kyr, specially for flux- 
cliangc thresholds larger than 50 % (see Figs. 8 and 9). This 
is also suggested from the analysis of the high temporal- 
resolution data in the next section. Although flux increments 
are more likely than decrements for any given threshold and 
time lag, a novel result of these simulations is that negative 
flux changes do happen, in contrast with the simple expec- 
tation of ever growing H II regions. 



3.5 Short-Term Variation Probabilities 

We re-ran four time intervals in each of Run A and Run 
B, spanning a few hundred years each, and producing data 
dumps at each simulation step (~ 10 yr). These time inter- 
vals were selected to contain a pair of negative/positive flux 
changes to investigate the correlation of flux variations with 
physical changes in the H II region, like size and density 
(next section). Therefore, they may not be representative of 
the entire sinmlation, but since it is not feasible to re- run the 
entire simulations producing data dumps at the highest tem- 
poral resolution, we use these data to constrain the expected 
flux variations in observable timescales. We argue that the 
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Figure 10. Probabilities for flux increments larger than 10 (top 
panel), 50 (middle) , and 90 % (bottom) as a function of time 
difference for the sample intervals at high time resolution. The 
error bars indicate the la statistical uncertainty from the number 
of counts in each bin 20-yr wide. 



1.0 



0.8 



a 0.6 



0.2 



0.0 



1.0 



0.8 



RunB 

Probability of flux increment > 10 % 



100 200 300 

Time Difference [yr] 



400 



500 



I" 0.6 



0.4 



0.2 



0.0 



RunB 

Probability of flux increment > 50 % 



******** t i i J 



100 200 300 

Time Difference [yr] 



400 



500 




200 300 
Time Difference [yr] 
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panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the sample intervals at high time resolution, for Run 
B. The error bars indicate the Icr statistical uncertainty from the 
number of counts in each bin 20-yr wide. 



close matcli at scales of At ~ 500 yr between the proba- 
bilities obtained from the low temporal-resolution (previous 
section) and the high temporal-resolution data (this section) 
indicates that the results presented here are meaningful. 

Figures 10 and 11 show, for Run A and Run B respec- 
tively, the probabilities for flux increments larger than the 
specified threshold as a function of time lag. These figures 
are the analogs of Figs. 6 and 7 for the high temporal- 
resolution data. The slight probability decrements after 
At ~ 300 yr, especially at low thresholds, arc an arti- 
fact caused by the fact that the data sets include a neg- 
ative/positive flux-change pair. Still, the probabilities at 
At = 490 yr match within 20 % to 80 % with the prob- 



abilities at At = 500 yr from the low-temporal resolution 
data. 

On observable timescales, At = to 40 yr, a small 
but non-zero fraction of H II regions is expected to have 
detectable flux increments. For two observations separated 
by 10 yr. Run A gives a prediction of 16.7 ± 2.9 % of H 
II regions having flux increments larger than 10 %, 6.8 ± 
1.9 % with flux increments larger than 50 %, and 4.7 ± 
1.6 % with increments larger than 90 %. The more realistic 
Run B predicts a smaller fraction of variable H II regions: 
6.9 ± 1.6 %, 0.3 ± 0.3 %, and % of them are expected to 
have flux increments larger than 10 %, 50 %, and 90 % over 
a time interval of 10 yr, respectively. 
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Figure 12. Probabilities for flux decrements larger than 10 (top 
panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the sample intervals at high time resolution. The 
error bars indicate the la statistical uncertainty from the number 
of counts in each bin 20-yr wide. 
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Figure 13. Probabilities for flux decrements larger than 10 (top 
panel), 50 (middle), and 90 % (bottom) as a function of time 
difference for the sample intervals at high time resolution, for Run 
B. The error bars indicate the Icr statistical uncertainty from the 
number of counts in each bin 20-yr wide. 



Figures 12 and 13 show the probabilities for flux decre- 
ments larger than a given threshold obtained with the 
high temporal-resolution data. The probabilities obtained at 
At = 490 yr from the high temporal-resolution data roughly 
match with those at At = 500 yr obtained from the low 
temporal-resolution data, within a factor of 1 to 3. 

Negative variations should also be detectable in a non- 
negligible fraction of H II regions. Run A predicts that 
5.7 ± 1.7 %, 2.6 ± 1.2 %, and 1.6 ± 0.1 % of H II regions 
should present flux decrements larger than 10 %, 50 %, and 
90 % respectively, when observed in two epochs separated 
by 10 yr. Run B predicts a smaller fraction of H II regions 
with negative flux variations: 3.3 ± 1.1 %, 1.5 ± 0.7 %, and 



% for thresholds at 10 %, 50 %, and 90 % respectively. 
We emphasize that negative variations arc not expected in 
classical models of monotonic H II region growth, while they 
are a natural outcome of the model presented here. More- 
over, these variations should be detectable with current tele- 
scopes. 

3.6 Variations in other properties of the HII 

regions 

Because Run B is more realistic in the treatment of 
fragmentation (Section 2.1, see Paper I for details), we use 
the time intervals of Run B with data at high temporal reso- 
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lution to investigate the correlations of sudden flux changes 
with other properties of the H II regions. 

Let the size of the H II region of interest be defined 
as Z/Hii = 2{A/'kY^'^ , where A is the area in the synthetic 
image where the H II region is brighter than three times 
the rms noise. Lets also consider the average density of the 
H II region pnii, and the rate at which neutral gas is con- 
verted into ionized gas A4"_>hii- Figure 14 compares the time 
evolution of these three quantities. 

Figure 15 further compares the time evolution of the 
2-cm fiux S'2cm, the total ionized mass Afnii in the H II 
region, and the denser ionized mass within the same volume 
Mhii, dense . The density threshold to define this dense gas 
is phii > 10~^^ g cm~^, the typical peak density reached 
when the H II region gets rapidly denser immediately after 
the fiickering events (see Fig. 14). 

The fiux-size correlation, as well as the size-density an- 
ticorrelation are a consequence of the relatively large optical 
depths of the H II regions (see Section 3.1). The "quench- 
ing" events, when an H II region has a large, sudden drop 
in flux and size, are coincident with large increments in the 
ionized density. The H II region rapidly reaches a state close 
to ionization equilibrium at its new size after the quenching 
instability. In Fig. 14 it is seen that at the moment of the 
quenching, the ionization rate has a sharp decrement imme- 
diately followed by an even faster increment that marks the 
initial re-growth of the H II region. Shortly afterwards, the 
ionization rate stabilizes again and the H II region grows hy- 
drodynamically, gradually becoming larger and less dense. 

The ionized mass of the H II regions follows the flux 
and size closely (Fig. 15), since this is the gas responsible 
for the free-free emission. Ifowever, if only the denser gas is 
taken into account, the amount of this dense ionized gas is 
particularly high in the initial re-growth of the H II region 
immediately after the quenching event (Fig. 15). Therefore 
the rapid quenching events are marked by the presence of 
denser gas around the accreting massive protostar. 



4 DISCUSSION 

4.1 A new view of early H II region evolution 

Until recently, H II regions have been modeled as 
(often spherical) bubbles of ionized gas freely expanding 
into a quiescent medium. However, this paradigm fails 
to explain observations of some UC and HC H II re- 
gions (see Section 1). There is enough evidence to assert 
that massive stars form in clusters by accretion of gas 
from thei r complex environment (reviews have been pre- 
sented by Garav fc Lizanol 19991: Mac Low fc Klessen|[200i : 
iBeuther et al. I l2007l : IZinnecker fc Yorkd |2007| ). Centrally 
peaked, anisotropic density gradients are expected at the 
moment an accreting massive star starts to ionize its envi- 
ronment. Therefore, the earliest H II regions should not be 
expected to be the aforementioned bubbles, but more com- 
plex systems where ionized gas that is outflowing, rotating, 
or even infalling may coexist. The fe asibility of this scenario 
has been shown by analvtic models |Hollenbach et al.l 1 19941 : 
lKetoll2002l . 12001 l2007l : iLugo et al.l l 2004^ and recent numer- 
ical simulations (Papers I, II, and III). In these numerical 
models, the inner part of the accretion flow is mostly ionized. 



while the outer part is mostly neutral. The neutral accretion 
fiow continuously tries to feed the central stars. The inter- 
action of this infalling neutral gas with the ionized region 
has a remarkable observational effect: the fiickering of the 
free- free emission from the H II region. 

We stress that the fiickering is not a gasdynamical ef- 
fect. Instead, it is a non-local result of the shielding of the 
ionizing source by its own accretion fiow. Hence, variations 
in the ionization state of the gas inside the H II region are 
not limited by the speed of sound but can happen on the 
much shorter recombination timescale of the ionized gas, 
rendering direct observation of this fiickering effect feasible. 



4.2 Observational Signatures 

Time-variation effects are expected in the ionized gas 
but not in the molecular gas, since small clumps of molecu- 
lar gas that become ionized and/or recombine have a much 
larger effect on the emission of the < 1 Mq of ionized gas 
than on the tens to hundreds of Mq of molecular gas in the 
accretion flow. One observati onal example is the HC H I I 
region G20.08N A, for which ICalvan-Madrid et al.1 (I2009D 
reports an ionized mass of 0.05 Mq and a mass of warm 
molecular gas in the inner 0.1 pc of 35 to 95 Mq (another 
well studied example is GlO.6-0.4, with more than 1000 Mq 
in th e pc-scale molecular flow surrounding the H II reg ions, 
e.g., iKetol Il990l : iKeto fc WoodI l2006l : iLiu et al.1 l201li . and 
references therein). The H II regions in the simulations also 
have masses of ionized gas in the range 10"'^ to 10~^ Mq 
while the stars are still accreting. 

In sections 3.4 and 3.5 we have attempted to quantify 
the expected flux variations in the radio continuum of H II 
regions for massive stars forming in isolation (Run A) and, 
more realistically, in clusters (Run B) . The accretion flow in 
Run A is stronger than in Run B (actually the star in Run 
A never stops accreting, see Papers I and HI), which leads 
to a brighter and more variable H II region. 

For a given run, flux variation threshold, and time dif- 
ference, positive variations are more likely to happen than 
negative ones, i.e., there is a constant struggle between the 
H II region trying to expand and the surrounding neutral gas 
trying to conflne it, with a statistical bias toward expansion. 

Monitoring H II regions for thousands of years is not 
possible, but a few observations of rapid flux changes over 
timescales of ^ 10 y r have been presented in the literature 
IIAcord et al. ' '1998*; 'Franco-Hernandez fc Rodriguez) |2004| : 
Ivan der Ta k ot al. 2005; Galvan-Madrid ct al. 2008). Com- 
paring multi-epoch images made with radio-interferometers 
can be challenging: even if the absolute flux scale of the 
standard quasars is known to better than 2 %, slight dif- 
ferences in observational parameters between observations 
(mainly the Fourier-space sampling, see jPer ley et al. igS^) 
make questionable any observed change smaller than 10 %. 
We have therefore measured the variation probabilities in 
the simulations at thresholds starting at 10 %. Considering 
this detection limit for the more realistic case of Run B, we 
have estimated that about 7 % of UC and HC H II regions 
observed at two epochs separated by about 10 years should 
have detectable flux increments, and that about 3 % should 
have detectable decrements. In total, ~ 10 % of H II regions 
should have detectable flux variations in a period of 10 years. 
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Figure 14. Variation at high temporal resolution in Run B of the 
scale length Lhii (filled black circles) of the H II region around 
the most massive sink particle, its average density Phii (filled 
black squares), and its rate of ionization of neutral gas M_>hii 
(filled black diamonds). 
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Figure 15. Variation at high temporal resolution in Run B of 
the 2-cm flux /S'2cm (filled black circles) of the H II region around 
the most massive sink particle, its ionized mass Mhii (filled black 
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Dedicated observations of as many sources as possible are 
now needed to test this model. 

Our long timescale data can only be constrained by ob- 
servational surveys, not by time monitoring. In Section 3.3 
we have shown that the radio luminosities (i.e., distance- 
corrected flux) in the long-term evolution of the simulated 
H II regions are consistent with major surveys, except for 
the most luminous H II regions, in which likely accretion 
has stopped and which therefore do not correspond with 
our data. The hypothesis that a considerable fraction of ob- 
served UC and HC H II regions may harbor stars that are 
still accreting material still needs more convincing evidence 
in addition to matching the model here presented. For most 
cases the dynamics of the surrounding molecular gas and of 
the ionized gas have not been studied at high angular reso- 
lution, and such studies in many sources are key to test this 
idea. From available observations, almost all of the massive 
star formation regions with signatures of active accretion 
and in which the mass of the protostar is estimated from 
dynamics to be > 20 Mq have a relatively bright H II 
regio n (with at least ~ 100 mjy at short cm wavelenghts, 
e.g.. iBeltran e t al. 2007i: iGalvAn^M adrid ct al. 2009). To our 
knowledge, the only clear exception is the recent report by 
IZapata et all (|2009( ) of an accreting protostar in W51 N with 
an estimated mass of M* ~ 60 Mq and only 17 mJy at 7 
mm. This object can be understood in the context of our 
simulations if it is in a quenched, faint state as currently 
observed. 



4.3 Caveats and limitations 

As mentioned in Papers I, II, and III, the simulations 
here presented do not include the effects of stellar winds and 
magnetically-driven jets originating from within 100 AU, 
which would produce outflows that are more powerful than 
the purely pressure-driven outflows that appear in Runs A 
and B (Paper I). 

The inclusion of stellar winds and jets may affect 
the results p resen ted in this paper only quantitatively. 
IPeters et al.l (|201ll ') have shown that magnetically driven 
outflows from radii beyond 100 AU do not stop accretion and 
even channel more material to the centr a l mos t massive pro- 
tostars. The simulations of I Wang et al. I (I2OIO') also indicate 
that coUimated outflows may be an important regulator of 
star formation by slowing the accretion rate but without im- 
peding accretion. Observationally, molecular outflows tend 
to be less coUimated for the more massive O-type protostars 
capab le of producing H II regions than for B-type proto- 
stars ijArce et al.ll2007l ). Regarding the radio-continuum, it 
is unknown if the free-free emission from the photoionized 
H II regions produced by O-type protostars can coexist with 
the free-free emission from (partially) ionized, magnetically 
driven jets. Before the appearance of an H II region, these 
jets are detected in protostars less massive than ~ 15 Mq 
(e.g.. ICarrasco-Gonzalez et al.ll2o"lol ). and even though their 
radio emission also appears to be varia ble (due to motion s 
and interactions with the medium, e.g.. lCuriel et al]|2006l ). 
their typical centimeter flux is ~ 1 mJy, one to two orders 
of magnitude fainter than the typical flux of UC and HC H 
II regions (except mayb e for the y oungest gravitationally- 
trapped H II regions, see lKetoll2003l ). Therefore, the relative 



effect of any variation in a hypothetical radio jet should be 
small compared with the variations in the H II region flux. 

A further limitation of this study is that accretion onto 
the protostars is not well resolved, since the maximum cell 
resolution (98 AU) corresponds to a scale of the order of the 
inner accretion disk (see also Paper I). 



5 CONCLUSIONS 

We performed an analysis of the radio-continuum 
variability in H II regions that appear in the radiation- 
hydrodynamic simulations of massive-star formation pre- 
sented in Paper I. The ultimate fate of ultracompact and hy- 
percompact H II regions is to expand, but during their evo- 
lution they flicker due to the complex interplay of the inner 
ionized gas and the outer neutral gas. The radio-luminosities 
of the H II regions formed by the accreting protostars in our 
simulations are in agreement with those of observational ra- 
dio surveys, except for the most luminous of the observed 
H II regions. We show that H II regions are highly variable 
in all timescales from 10 to 10* yr, and estimate that at 
least 10 % of observed ultracompact and hypercompact H 
II regions should exhibit flux variations larger than 10 % for 
time intervals longer than about 10 yr. 
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